Accepted for publication in The Astrophysical Journal 

Preprint typeset using LAT^X style cmulateapj v. 10/09/06 



COSMIC REIONIZATION AND THE 21-CM SIGNAL: COMPARISON BETWEEN AN ANALYTICAL MODEL 

AND A SIMULATION 

Mario G. Santos 1 , Alexandre Amblard 2 , Jonathan Pritchard 3,4 , Hy Trac 4,5 , Renyue Cen 5 , Asantha Cooray 2 

1 CENTRA, Institute Superior Tecnico, Lisboa 1049-001, Portugal 
2 Center for Cosmology, Department of Physics and Astronomy, 4186 Frederick Reines Hall, University of California, Irvine, CA 92697 

California Institute of Technology, Mail Code 130-33, Pasadena, CA 91125 
4 Harvard-Smithsonian Center for Astrophysics, MS-51, 60 Garden St, Cambridge, MA 02138, USA 
°Dcpartmcnt of Astrophysical Sciences, Princeton University, Princeton, NJ 08544 
Draft version August 28, 2008 

ABSTRACT 

We measure several properties of the reionization process and the corresponding low- frequency 21- 
cm signal associated with the neutral hydrogen distribution, using a large volume, high resolution 
simulation of cosmic reionization. The brightness temperature of the 21-cm signal is derived by 
post-processing this numerical simulation with a semi-analytical prescription. Our study extends 
to high redshifts (z ~ 25) where, in addition to collisional coupling, our post-processed simulations 
take into account the inhomogencitics in the heating of the neutral gas by X-rays and the effect 
of an inhomogeneous Lya radiation field. Unlike the well-studied case where spin temperature is 
assumed to be significantly greater than the temperature of the cosmic microwave background due to 
uniform heating of the gas by X-rays, spatial fluctuations in both the Lya radiation field and X-ray 
intensity impact predictions related to the brightness temperature at z > 10, during the early stages of 
reionization and gas heating. The statistics of the 21-cm signal from our simulation are then compared 
to existing analytical models in the literature and we find that these analytical models provide a 
reasonably accurate description of the 21-cm power spectrum at z < 10. Such an agreement is useful 
since analytical models arc better suited to quickly explore the full astrophysical and cosmological 
parameter space relevant for future 21-cm surveys. We find, nevertheless, non-negligible differences 
that can be attributed to differences in the inhomogeneous X-ray heating and Lya coupling at z > 10 
and, with upcoming interfcrometric data, these differences in return can provide a way to better 
understand the astrophysical processes during reionization. 

Subject headings: cosmology: theory — large scale structure of Universe — diffuse radiation 



1. INTRODUCTION 

One of the key challenges faced today in cosmology is 
to understand in detail how the density distribution of 
both dark matter and baryons in the Universe evolved 
from a relatively smooth initial state at early times into 
the non-linear structures we observe today. This non- 
linear structure formation is directly coupled to the for- 
mation of galaxies first and later, galaxy clusters. The 
epoch of reionization (EoR) is a crucial stage in the his- 
tory of galaxy and structure formation, signaling the 
birth of the first luminous objects as structures first 
evolved beyond the well understood linear regime. Al- 
though the process by which the intcrgalactic medium 
(IGM) became ionized is quite complex, the current view 
is that, when the first proto-galaxies and quasars form, 
they ionize the surrounding gas creating the HII "bub- 
bles". These regions continue to grow and overlap, so 
that eve ntually all of the neutral gas in the IGM b ecomes 
ionized (jBarkana fc Loebll2001t IFan et alJl2006af ). 

Current primary constraints on the epoch of reion- 
ization come from two main sources: the Wilkinson 
Microwave Anisotropy Probe (WMAP) determination 
of the optical depth to recombination through a late- 
time signature at large angular scales in the cosmic 
microwave background (CMB) pola rization spectrum 
(jSpergel et alJl2006t IZaldar riagd[l997l ) and the Lya for- 
est absorption spectra t owards high redshift quasars (e.g. 
IFan et al.ll200l l2006b[ ). This latter Gunn-Peterson ef- 



fect Jg unn fc PetersonI |1965[ ) is pr esent towards sight 
lines to quasars out t o z ~ 6.5 (jBecker et al.l 120011 : 
ICen fc Mc Donald 2002) showing that reionization should 
be ending by this time (it can also mean that there is a 
transition in the Gunn-Peterson optical depth from ab- 
sorption spectra out to z ~ 4 and those out to higher 
redshifts). However, we note that a small neutral frac- 
tion is enough to completely absorb the Lya quasar flux, 
so these observations themselves cannot be used to prop- 
erly establish th e reionization history of the Universe 
(|Lidz et alJl2006h . 

In terms of t he WMAP data , the large angular 
scale polarization (jPage et al.ll2006HDunklev et al.ll2008ft 
yields a Thomson optical depth of r = 0.084 ± 0.016, so 
that reionization should happen at z ~ 10.8 ±1.4 in 
the favored ACDM cosmological model, i f we assume in- 
stantaneous reionization of the Universe (jKomatsu et alJ 
120081 : ISpergel et al.|[2006l ) . The reionization process need 
not be instantaneous a nd if it is less abrupt (see e.g. 
lHaiman fc Hol der 2003) then the Universe may have be- 
gun to partly reionize at an even earlier epoch. Lim- 
ited by these two constraints, we know that the reioniza- 
tion process should have lasted for at least 500 million 
years, although there is very little observational evidence 
on how this event actually occurred, allowin g for vari- 
ous possibilities includ ing double reionization (|Cenll2003l : 
IWvithe fc Loeb1l200lh . 

While more precise CMB polarization measure- 
ments than with WMAP alone at large angular 



scales can p rovide more information on the reioniza- 
tion history (iHolder et al.l l2003l: iKaplinghat et al.ll2003l : 
iMortonson fc Hul l2008f ). it is now generally accepted 
that detailed information, including the exact his- 
tory of the reionization process, will become available 
with 21-cm signal from the neutral hy drogen distri- 
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Given the line emission, with frequency selection for 
observations, the 21-cm data pro vide a tomographic 
view of the reionizati on process (|Santos et all 120051 : 
iFurlanetto et al.ll2004bO as well as a probe of the dark 
ages where no luminous sour ces are present after recom - 
bination at a redshift of 1100 (|Loeb fc Zaldarriagall2004T ) . 

We note that small angular scale CMB anisotropics 
also capture some information related to reionization, 
especially the inhom ogeneous or patchy nature of the 
reionization process (ISantos et al. 2003; Aghani m et all 



Il996t iKnox et all I1998T) and through effects s u ch as the 
Ostri ker-Vishniac ( Ostriker fc Vishniad 119861 ; IVishniad 
119871 ) effect. The 21-cm background and CMB provide 
complimentary information related to reionization since 
the former is related to the neutral hydrogen dist ribution 
while the latter is due to the electron content (jCooravl 
2004). Unfortunately, at such small angular scales, the 
CMB anisotropy spectrum is rich with a variety of ef- 
fects contributing to the overall signal, including galaxy 
clusters and gravitational lensing. Therefore, the focus 
is mainly on 21-cm observations, but additional infor- 
mation from CMB, such as with a high-resolutio n ver- 
sion o f the EPIC concept mission for the CMBpol (|Boc1 
120081 ) may help extract some properties of the reioniza- 
tion physics. 

Motivated by the existing observational constraints 
and the possibilities to study reionization through the 
neutral hydrogen content with the proposed 21-cm ex- 
periments such as the Square Kilometer Array (SKA 1 ), 
the Low Frequency Array (LOFAR 2 ) and the Milcura 
Wide-field Array (MWA 3 ), a great deal of effort has 
been made recently to unders tand the 21-cm signal 
and its information content (see IFurlanetto et alll2006bl 
for a review). In parallel with developments in the 
experimental front, our theoretical understanding of 
reionization has also improved both through numeri- 
cal simulations and analytical models. Numerical sim- 
ulations provide a detailed description of related astro- 
physics at these rcdshifts from first principles by directly 
solving the non-linear physics of gravita tional co l lapse , 
hydrodynamics, and radiative transferdGnedinl 120001: 
Razoumov et alll2002l ; ISokasian et all l l200aiCiardiet"al 
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2001: iMcQuinn et alll2007l i However, proper sampling 
of the epoch of reionization re quires simulations with 
large volumes (~ 100M pc/h) 3 ((Barkana fc Loebl 12004 ; 
IFurlanetto et all l2004d ) and with adequate mass reso- 
lution to resolve halos where first-light sources are ex- 
pected to form (M ~ 10 6 M Q /h). The large volume and 
the large particle number makes such simulations compu- 

1 http://www.skatelescope.org 

2 http://www.lofar.org 

3 http:/ /www. haystack. mit.edu/arrays/MWA 



tationally expensive, especially in the context of hydro- 
dynamical calculations related to the gas physics. The 
usual approach is to use high resolution, but small vol- 
ume simulations or large volume but low mass resolution 
simulations. In this paper we make use of a simulation 
that directly resolves halos below 10 s M in a box of 100 
Mpc/h, considered essential to properl y take into accoun t 
the bubble growth during reionization (jShin et al.ll2007l ). 

Due to the challenging computational requirements of 
numerical models, progress on the modeling front has 
come mostly from analytical descriptions on the volume 
filling factor, size distribution of ionized regions, as well 
as the pow er spectrum of the ioniz e d fract i on and den- 
sity fields (IFurlanetto et~aTI l2004d . l2006at ISethil 120051 ; 
Barkana 2007). These analytical descriptions have been 
quite useful to understand the p ossible contributions to 
the 21-cm signal at high redsh ifts ([Barkana fc Loeb1l2005l ; 
iPritchard fc Furlanettoll2007f ) or under certain simplified 
conditions such as the case where spin temperature of 
neutral gas is significantly higher than that of the CMB. 
The analytical approach is also crucial to explore the 
extent to which the full parameter space of the 21-cm 
signal and associated cosmology can be established with 
data from future surveys planned with MWA, LOFAR , 
and SKA (ISantos fc Cooravl 120061 ; IMcQuinn etaTl 120061 ; 
iMao et al.ll200an . 

An intermediate approach, based on semi- 
analytical models combined wit h semi-numerical 
models, has also been de veloped (jZahn et all 120071 : 
iMesinger fc Furlanettoll2~007f ). It relies on the generation 
of realizations of halo distributions directly from the 
linear density field and implementing the corresponding 
ionization map using criteria similar to the analytical 
models. These allow to preserve the spatial information 
of the reionization process as provided by simulations, 
while achieving a much larger dynamic range than 
provided by radiative transfer codes. 

In future, once data become available with the first- 
generation low- frequency radio interferometers, it will 
be useful to have fast techniques to extract the parame- 
ters from the measurements, such as the power spectrum 
21-cm brightness temperature. While detailed numeri- 
cal simulations to scmi-numcrical models can be consid- 
ered for this purpose, it is unlikely such simulations can 
be carried out for all variations in parameters of inter- 
est, which includes both astrophysical and cosmological 
quantities. In this sense, it is more useful to make use of 
analytical methods supplemented by well-motivated fit- 
ting functions for quantities such as the power spectrum 
of ionization fraction during reionization in terms of the 
power spectrum of density perturbations that depends 
on cosmological parameters. This is the approach taken 
in predi ctions related to 21-cm cosmological info rmation 
content (|Santos fc Cooravll2006l : lMao et all2008h . but we 
still need to improve our approximations in such an ap- 
proach by continuing to study first-principle numerical 
models of reionization and 21-cm physics to test assump- 
tions on the existing analytical models. 

Due to time constraints associated with numerical 
models, it is very likely that analytical models are the 
preferred choice for intensive astrophysical and cosmo- 
logical parameter studies from the 21-cm signal observed 
with low-frequency radio interferometers. In the case of 
cosmological parameter estimation with CMB or from 
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galaxy clustering, the numerical computations provide 
solutions to an analytical derivation of either the CMB 
signal or the dark matter clustering power spectrum. 
Unfortunately, for the 21-cm anisotropics, such an ap- 
proach is likely to be complicated, especially during 
rcionization though at very high rcdshifts (~ 50), where 
physics is simple, a quick numerical calculation of the 
physics involved can be c arried out from first-principles 
([Lewis fe Challinori[200l . 

In this paper we determine several properties of rcion- 
ization using a state of the art large volume and high- 
resolution simulation of cosmic rcionization based on 
a photon-advection radiative transfer scheme combined 
with a dark matter N-body simulation with recipes for 
baryons and star formation (|Shin et al.ll2007f ) . While the 
simulation itself is dark matter, the post-processing al- 
lows us to convert the star-formation rate predicted in 
the simulation to the spin temperature of the gas in 
the simulation. Furthermore, we determine the 21-cm 
brightness temperature up to z ~ 25, by post-processing 
the simulation output with a semi-analytical prescription 
for the X-ray heating of the gas, the Lya coupling, and 
the collisional coupling, and by taking fully into account 
the spatial fluctuations in these quantities. The existing 
simulations of the reionization and predictions related 
to 21-cm signal from such simulations generally ignore 
the spatial inhomogeneities associated with X-ray heat- 
ing and Lya radiation field, though anisotropies of the 
21-cm brightness temperature is expected to be sourced 
by such inhomogeneities at high redshifts. We compare 
results based on simulations wit h estimates from a fast 
analytical model of reionization ( Furlanetto et al.ll2004d : 
IMcQuinn et al.l 120051 : lUdzet al.ll2007D~ 

Throughout the paper, we make use of the follow- 
ing cosmological parameters: il m = 0.26, Qa = 0.74, 
fib = 0.044, h = 0.72 and n s = 0.96, ct 8 = 0.77 based on 
the results from W MAP, SDSS, BAO, SN and HST (sec 
iSpergel et al.l [20061 and references therein). The optical 
depth is r « 0.09, consistent with the WMAP-5 result. 
The paper is organized as follows: In the next Section, 
we outline details related to the rcionization process and 
compare results from the simulation and the analytical 
calculation. We then proceed to describe how to cal- 
culate the corresponding 21-cm signal in Section 3 with 
details of the simulation in Section 4. Again we show a 
comparison of the results from simulation to analytical 
models (Section 5). We conclude with a summary of our 
results in Section 6. 

2. COSMIC REIONIZATION 

2.1. Numerical Simulation 

In this paper, we make use of one of the largest simu- 
latio ns of cosmic reion ization that has been co mpleted to 
date ([Shin et alj|2007t ). We refer the reader to lShin et all 
(|2007[ ) for details related to the hybrid code that con- 
tains a N-body algorithm for dark matter, prescriptions 
for baryons and star formation, and a radiative transfer 
(RT) algorithm for ionizing photons. We provide a basic 
summary of the simulation parameters here as necessary 
for this study. 

The hybrid simulation involves a high resolution N- 
body calculation of 2880 2 dark matter particles in a 
L = 100 Mpc/h box. With a particle mass resolution 
of 3.02 x 10 6 Mq/H, halos can be reliably resolved down 



to masses of ~ 10 8 M Q //i, accounting for the majority 
of photo- ionizing sources. The simulation distinguishes 
between the first generation, Population III (Pop-Ill) 
stars and the second generation, Population II (Pop- 
II) stars by following the ch emical enrichment of the 
ISM and IGM as described in iTrac fc Ce"nl (f2007h . The 
input UV spectrum is divided in three energy ranges 
13.61 eV< hv < 24.59 eV, 24.59 eV< hv < 54.42 eV 
and hv > 54.42 eV, with Pop-II stars with a Salpeter 
IMF providing 1100, 3830 and 270 ionizing photons per 
baryon of star formation respectively. For Pop-Ill stars 
with a top-heavy IMF, the correspond i ng nu mbers are 
15000, 51500, and 3500 (|Schaeredl200l . [20031) . The ra- 
diative transfer of ionizing photons is calculated simul- 
taneously as dark matter evolves with the N-body code 
and with star formation and baryon physics evolving ac- 
cording to recipes each step of the way. In this way, 
our simulation differs from other descriptions in the lit- 
erature where radiative transfer and baryon physics are 
obtained by post-processing a completed N- body run. 

Not e that we do not use the halo model of lTrac fc Cenl 
(|2007| ) for prescribing baryons and star formation. In- 
stead, an alternative approach is taken where we calcu- 
late the local matter density p and velocity dispersion a v 
for each particle. The baryons are assumed to trace the 
dark matter distribution on all scales and we obtain the 
local baryon density p^ = p(Q(,/D, m ) and gas tempera- 
ture T = po~ v 2/(3k). Star formation is only allowed to 
occur in particles with densities p > 100p C rit(z) and tem- 
peratures T > 10 K, thus restricting star formation to 
regions within the virial radius of larger halos and these 
halos are fully resolved given the low mass resolution of 
our simulation. 

The radiative transfer of ionizing radiation uses a 
photon-advection scheme and was run simultaneously 
with the N-body calculations using a RT grid with 360 3 
cells. However, the ionization and recombination calcula- 
tions were done for each particle individually rather than 
on the grid to preserve small-scale information down to 
scales of several comoving Kpc/ft.. For post-processing, 
the dark matter, baryons, and radiation are collected on 
a grid with 720 3 cells and the data are saved every 10 
million years from z = 25 down to z = 5. 

We start our analysis of the reionization process by ex- 
amining a sequence of cuts through the simulation box 
of the fraction of free electrons x\ (the ionization frac- 
tion), shown in Figure [T] In this simulation box, reion- 
ization begins roughly at a redshift of ~ 18. The ionized 
bubbles have complex shapes and cannot be simply de- 
scribed with spherical models for the 3-dimcnsional HII 
regions surrounding UV sources. Complete overlap of 
the ionization patches occurs at the redshift of z ~ 6. 
Note that although the ionization fraction can have a 
range of values between and 1, most of the volume in 
the simulation is almost completely ionized (very high 
Xi) or completely neutral, thus favoring the current view 
of reionization based on the percolation of large ionized 
bubbles. Figure shows the evolution of the average ion- 
ization fraction compared with the values if we assume 
that gas is completely ionized inside bubbles (xi = 1) 
and completely neutral outside (xi = 0). Bubbles are 
defined by the threshold X{ > 0.5. 

2.2. Analytical models 
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Fig. 1. — Maps of the ionization fraction from the simulation at rcdshifts z = 20.60,15.24,10.00,7.40, corresponding to Xi = 
0.0002, 0.03, 0.35, 0.84. Note how there is a clear separation between the highly ionized regions (in red) and the mostly neutral IGM 
(black). 

Our analytic model for reionization, which we com- 
pare with result s from our numer i cal sim ulation, follows 
the approach of iFurlanetto et all j2004c). The mass of 
the ionized gas is linked to the mass in galaxies by the 
ansatz, m ion = C m gah where £ is an ionizing efficiency. 
A spherical region of gas of mass m is considered ion- 
ized if it contains sufficient sources to self ionize, i.e. 
/coil > C • I n the excursion set formalism this crite- 
ria is well described by a mass dependent linear barrier 
B(m,z) = Bq + B\o 2 {m), where a 2 (m) is the variance 
of the density fluctuations on the scale m. With this we 
can calculate the mass function of bubbles (the comoving 
number density of HII regions with masses in the range 
m ± dm/2): 



dn(m) 12 p 

dm V 7r m 2 



dlogcr 



dlogm 



a(m) 



exp 



B 2 (m,z) 



2a 2 (m) 



dm . 
(1) 

where p is the mean mass density of the Universe. 
Note that we renormalizc the resulting mass function 



to enforce the requirement Q = C/coilj where Q is the 
filling fraction of bubbles. Next we must determine the 
appropriate value for the ionization efficiency £. We al- 
low £ to vary with redshift and require that Xi = C/coii 
at all redshifts. Here Xi is determined from the simula- 
tions so that we bypass the need for a source prescrip- 
tion when defining Q and we assume a Press-Schechter 
mass function when determining / co ii- In principle, 
using the Sheth-Tormen mass function gives a weaker 
redshift dependence for £, but we use Press-Schechter 
for greater consistency with the reionization model of 
IFurlanetto et all (|2004ch . 

To calculate fluctuations in the 21-cm brightness tem- 
perature, we first need to calculate the correlation func- 
tions in the ionization fraction £ XiXi , density £gg, and the 
cross-correlation of these two quantities £ Xi 5, where 



U = ((a-a)(b-b)) 



(2) 



and S = p/p — 1. We use th e halo model to calculate 
^ (jCoorav fc Shethl l2002t ). IFurlanetto eTaTI (|2004cf ) 
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Fig. 2. — The average ionization fraction as a function of redshift 
from the simulation (red solid line) and assuming complete ioniza- 
tion of the bubbles (blue dashed line), defined for x; > 0.5. The 
green dot-dashed line shows the ionized fraction inside the bubbles 
while the cyan dotted line shows the ionized fraction in the IGM 
(defined has Xi < 0.5). 

present an ad hoc model for the correlation functions £ xx 
and £ x g, designed to ensure that the correct limiting be- 
havior as xh — ► 0, 1 is obeyed. A fundamental problem 
with their approach is that bubbles arc assumed to be 
spherical at all times leadin g to problems describin g the 
overlap of bubbles prope rly. iMcQuinn et alj (12005T) later 
attempted to modify the lFurlanetto et al.l ( 2004cD model 
to forbid bubble overlap. However, since neither of these 
models correctly handles bub ble overlap, we choose t o 
use the original formulation of iFurlanetto et ail (|2004cD . 
However, we incorporat e the corrected calcula tion of the 
bubble bias, as noted by McQui rm et al.l (|2005f ). We note 
that a more physically motivated method based upon the 
two-step appro ximation has recently been developed by 
iBarkanal (|200l . but for purposes of the present discus- 
sion where we are investigating the extent to which a 
simple model for parameter estimates can be compared 
to numerical simulations and is found in general to pro- 
vide good agreement, such details are unnecessary given 
the availability of simulations. 

The 3-dimensional power spectra of our simulation 
were performed using the fast Fourier transforrn pack- 
age fftw-3.1.1 4 and is defined through: < a(k)b*(k') > = 
(27r) 3 <5 3 (fc — k')P a i,(k). We then binned our modes with 
Sk = ^(Mpc/Zi)^ 1 and computed the average power 
spectra in each bin. Throughout the paper we will plot 
the dimensionless power spectrum, A 2 = fc 3 P(fc)/27r 2 , 
which gives the contribution to the variance per logarith- 
mic interval in k. In order to test the analytical calcula- 
tion, we show in Figure [3] the ratio of the 3-dimcnsional 
power spectrum of the ionization fraction to the matter 
density, for the simulation and analytical models. 

4 |http: / /www.fftw.org| 



The features in the power spectrum of ionization frac- 
tion relative to density perturbations can be described 
as following. At scales much larger than the bubble size, 
the ionization fraction power spectrum is proportional to 
the matter density with an overall scaling which can be 
assigned to the bias factor of the ionized regions. From 
the simulation results shown in dotted lines this can be 
seen more easily at high redshifts where bubble sizes are 
small. At these redshifts (z = 20.60 and z = 15.24), 
there is an increase in the ionization power spectrum at 
small scales due to the Poisson behavior of the bubble 
distribution. This is analogous to the shot-noise compo- 
nent of the galaxy power spectrum even at low redshifts 
that dominate fluctuations at non- linear scales. 

As we move to lower redshifts and larger bubbles, Fig- 
ure [3] shows more clearly that the power spectrum of 
the ionization fraction P XiXi decreases on physical scales 
smaller than the typical bubble size, due to the smooth- 
ing effect of bubbles. Note however that the decrease is 
not abrupt since there is a distribution of bubble sizes 
at any redshift with varying ionization fractions. The 
analytical calculation we discussed so far (dashed line) 
seems to agree reasonably well with the simulation, al- 
though at larger ionization fractions (xi > 0.8), when the 
bubbles occupy most of the simulation volume, there are 
some differences. This difference is probably due to the 
finite size of the simulation box which effectively limits 
the bubble size and reduces power on large scales. Note 
however that even if we our simulations involved larger 
volumes, the power could continue to be smaller than 
in the analytical case since the simulation accounts for 
sclf-shiclding of dense regions that remain neutral. This 
causes the bubble gro wth to stall and in return limit the 
maximum bubble size (jFurlanetto fc Ohll2005f ). 

In Figure [31 in addition to the comparison between re- 
sults from our numerical simulation and a widely-used 
analytical model, we also compare with two simple fit- 
ting functions of the ionization fraction given the density 
field power spectrum. We do this comparison since when 
exploring t he full parameter space probed by 21-c m ex- 
periments (jSantos fc Cooravll2006HMao et al.ll2008l ). cal- 
culations make use of fitting formulae to describe the 
power spectrum of the ionization fraction without mak- 
ing use of detailed analytical models or numerical calcu- 
lations of bubble growth. Here, we compare our results 
to two fitting functions in the literature where the power 
spectrum of the ionized fraction P XiXi is given as 

-{kR x 



P, 
P, 



*«(*) 



C t 1 + 



Pss(k) 
ti {kR XiXi )-\- 

- 7 /2 



(3) 
(4) 



+ (fciW 2 ] "* Pss(k) 
where the first approximation is from | Santos fc Cooravl 
(|2006t ) and the second is from iMao et al.l (|2008[ ) and 
bxiXii &xiXi, 7, and R XiXi are free parameters that are 
varied to obtain a fit to the simulation results. In Ta- 
ble 1, we list the values that were obtained by comparing 
numerical simulations with the above form. 

While these two fitting functions have been used in the 
literature when making predictions related to how well 
cosmological parameters can be measured with 21-cm in- 
terferometers such as MWA and LOFAR, as can be seen 
from Figure [3J both these functions are not completely 
accurate descriptions of the power spectrum of ionized 
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Fig. 3. — Ratio of the power spectra of the ionization fraction, 5 Xi = x%/xi — 1 to the matter one, <5 = p/p — 1. The simulation 
corresponds to the dott ed (black) line while the dashed-green line gives the results of the analytical model followe d in the paper (based on 
IFurlan ctto et al. 2004c). Dot-dashed (b lue) lines a r e fitti ng curves using the two parameter model suggested in ISantos fc Cooravl ||2006D 
while solid (red) curves use the model in Mao ct al. (2008). Rcdshifts, as labeled from top-left to bottom-right are: z = 20.6, 15.2, 10.0, 7.4, 
corresponding to Sii = 0.0002, 0.03, 0.35, 0.84 respectively. 



Function Parameter z=20.6 z=15.2 z=10.0 z=7.4 

Eq.[3] R x%Xi (h^Mpc) 0~ 016 042 0.96 

b XiXi 27 17.7 4.0 1.3 

Eq.[5] R XzXt (h-^Mpc) -4.53 0.23 0.62 1.24 
a x%Xi 0.83 -1.12 0.93 1.48 
7 -1.77 3.72 1.93 1.99 
b XiXi 27.2 13.1 4.7 1.6 



TABLE 1 

Parameters of the ionization fraction power spectrum 

model described in equations [3] & \e\ fitted to our 
simulation power spectrum measured at 2= 20.6, 15.2, 10, 
and 7.4. Note that to get the "physical" R XiXi based on 
our definitions one should multiply the above values by a 

FACTOR OF (27t). 



fraction over all the rcdshift range we have studied with 
simulations. This is due to the fact that, for exam ple, 
the first approximation from lSantos fc Cooravl (|2006f ) as- 
sumes that bubbles can be described with a single size 
leading to a sharp cut-off at the wavenumber correspond- 
ing to the inverse radi us, while from simulatio ns and in 
the analytical model of lFurlanetto et al. ( 2004c]) the bub- 
bles have varying sizes. At z < 15, when bubbles have 
started to grow, the fitting formula of iMao et alj (|2008f ) 



provides an improved fit given the additional freedom 
provided by the parameters a and 7 and this descrip- 
tion is probably adequate enough for now when making 
predictions related to the extent to which cosmological 
parameters can be measured with a 21-cm experiment. 
On the other hand, when model fitting real data, it may 
be necessary to improve the estimate of P Xi xi (k) beyond 
simple fitting functions as the ones listed in equations (|3j) 
and (El) . In this respect , we no te that the analytical 
model of lFurlanetto et all (|2004cD provides a more accu- 
rate description of the results from our numerical sim- 
ulation and a more clear interpretation of the values 
obtained. If such a model can be further improved to 
quickly explore a large parameter space in a reasonable 
time, it may be useful to implement such a model in a 
numerical code for parameter estimates, such as based 
on the Markov Chain Monte-Carlo technique, instead of 
simply using fitting functions to estimate parameter val- 
ues from 21-cm interferometers. 

Figure 2J left panel, shows the cross power spectrum 
between the ionization fraction and the density field. The 
peak of this cross-power spectrum moves to larger scales 
as the redshift decreases since it is related to the typical 
size of bubbles during reionization. At smaller scales, 
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Fig. 4. — Left: Cross power spectra of the ionization fraction and matter density ((fc^)) from simulations. Solid lines give the 
simulation values while dashed lines are for the analytical model from Furlanetto et al. From top to bottom the curves correspond to 
z = 20.6,15.2,10.0,7.4. Right: The correlation coefficient, r = P Xi s/ ' *fPx~x~Ps5, with simulations (solid line) and the fitting function of 
IMao et al.l 120081) (dashed line) for the same redshifts as the left panel. 



there is some indication that the cross-correlation power 
spectrum becomes negative in the numerical simulation 
suggesting more of an anti-correlation than what is seen 
in the analytical model. We believe this is partly due 
to small neutral regions inside HII bubbles because the 
simulation takes into account the self-shielding of dense 
regions. Note however, that as reionization progresses 
and the radiation intensity become larger, the existence 
of small fingers from large bubbles protruding into the 
now small neutral regions will play an important part in 
this anti-correlation. 

Similar to equation ([5]), we can also write the cross 
power spectrum between Xi and 8 assuming a perfect 
correlation of the two fields such that 

/', , x /', , /',, . (5) 
and by making use of the fitting formulae for the ion- 
ization fracti on power spectrum from equat ion (|3"|). as 
was u sed by ISantos fc Cooravi (|2006l ). In IMao et aTl 
(|2008f ). this cross power spectrum is modeled as P Xi $ = 
bl zS exp[-a XtS (kR XiS ) - (kR XtS ) 2 ]Pss, and in Figured 
right panel we also show this case with a new set of pa- 
rameters (b Xi s,a Xi s, R Xi s). Here, we plot the correlation 
coefficient r = P Xi $/ \J P XiXi P$s and we compare the fit- 
ting function motivated bv lMao et al.l (|2008h to numeri- 
cal simulations. While there is a good overall agreement 
between the simulations and the analytical model, we 
find differences at z ~ 20 between the two descriptions. 
Such a disagreement, however, is not a concern for first- 
generation 21-cm observations since these instruments 
mostly concentrate on z ~ 6 to 9 during reionization. 

3. THE 21-CM SIGNAL: THEORY 

We present here the overall calculation of the 21-cm 
signal. Details specific to the simulation and the ana- 
lytical calculation considered in this paper will be given 
later in the appropriate sections. 



3.1. Brightness temperature 

One of the best ways to observe the reionization pro- 
cess is through the 21-cm brightness temperature, cor- 
responding to the change in the intensity of the CMB 
radiation due to absorption or emission when it travels 
through a patch of neutral hydrogen. It is given, at an 
observed frequency v in the direction n, by 



ST b (n,u) 



1 + z 



(6) 



where Tg is the temperature of the source (the spin tem- 
perature of the IGM), z is the redshift corresponding 
to the frequency of observation (1 + z = v<x\jv^ with 
i/ 2 i = 1420 MHz) and T 7 = 2.73(1 + z)K is the CMB 
temperature at redshift z. The o ptical depth , r, of this 
patch in the hyperfine transition (|Fieldl[l959] ) is given in 
the limit of ksT s » hf2i by 

3c 3 HA 10 rim ,„s 

T ~ l&ku^Tsil + z^dVr/drY [) 

where A±o is the spontaneous emission coefficient for 
the transition (2.85 x 10~ 15 s" 1 ), tt-hi is the neutral hy- 
drogen number density and dV r /dr is the gradient of 
the total radial velocity along the line of sight (with 
V r = V • n); on average dV r /dr = H(z)/(1 + z). In this 
paper we will neglect perturbations from the peculiar ve- 
locity of the gas. The neutral density can be expressed 
as n m = fmXp b /m p where f HI = pm/pti is the frac- 
tion of neutral hydrogen (mass weighted), X ps 0.76 is 
the hydrogen mass fraction, p\, is the baryon density and 
m p the proton mass. The 21-cm temperature is then: 

0.7, 



6T b (n,v) 



;23/ H i^ 
Pb 



1-*L 



n b h 2 

0.02 



0.15 



1 + z 
10 



"I 1/2 



mK (8) 
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In order to proceed, we will need a prescription to calcu- 
late the spin temperature of the gas. 

3.2. Spin temperature 

The spin temperature is coupled to the hydrogen gas 
temperature (Tk) through the spin- flip transition, which 
can be excited by collisions or by the absorption of Lya 
photons (Wouthusysen-Field effect) and we can write: 



1 



T s 



ytot 



where ytot = Ua + Vc is the sum of the radiative and colli- 
sional coupling parameters and we arc already assuming 
that the color temperature of the Lya radiation field at 
the Lya frequency is equal to Tk- When the coupling to 
the gas temperature is negligible (e.g. y tot ~ 0), Ts ~ T 7 
and there is no signal. On the other hand, for large ytot-, 
Ts simply follows Tk- 

Collisions can be important for decoupling the HI 21- 
cm spin tempera t ure fr om the CMB, especially at high 
rcdshifts iNusserl (|2005) and the coupling coefficient is 
given by 

4T* 



Vc = 



SAinT-f 



[ K ? H (T k )n H + K l H (T k )n e 



(10) 



where = fec/fcA 2i C m = 0.0628 K, Atf_^ is tabulated 
as a function of T k $M ison fc Dalgarnolll969l ; IZvgermanl 
l2005h . 

K i-o m taken from iFurlanetto fc Furlanettol 
(1200% and ri e is t he electron number density (see also 
iKuhlen et aT1l2006t ). For a more detailed analysis of the 

collisional coupling, see IHirata" fc Sigurdson] (120071) 

Th e Wouthvsen-Field effect ([Wouthuvse'r] 119521; iFieldl 
fl959h coupling is given by 



~17 



with 



167r 2 r,e 2 / a 
27A 10 T 7 r 



\T^m e c 
;5.552 x 10" 



(11) 



(12) 



(1 + *) 



Hz 



where f a = 0.4162 is the oscillator strength of the 
Lya transition and (1 + z) comes from T 7 . In equa- 
tion 111! S a is a correction factor of order unity, 
which describes the detailed structure of the pho- 
ton distr ibution in the neighborhood o f the Lya res- 
onance (IChen fc Miralda-Escudel 12004: IHirata! 120061 : 
IChuzhov fc Shapirdl2007l : IFurlanetto fc Pritchar J2006I ) . 
We make use of the ap proximation for S a outlined in 
IFurlanetto etal] (|2006bT ). The proper Lya photon inten- 
sity, Jo, (the spherical average of the number of photons 
hitting a gas element per unit proper area per unit time 
per unit frequency per steradian) is given by a sum over 
the hydrogen levels n, 



J a (x,z) 



(1 + ^) 2 

47T 



(13) 



n=2 



I — J dx' e a (x + x',L>' n ,z') , 

where frec{n) is the fraction of Lyman- n photon that cas- 
cade through Lya and e(x, v n , z) is the comoving photon 



emissivity, defined as the number of photons emitted at 
position x, redshift z and frequency v per comoving vol- 
ume, per proper time and frequency. Note that the red- 
shift z 1 in equation [T4l is such that x' = /* cH~ 1 dz". 
The absorption at level n at redshift z corresponds to an 
emitted frequency at z' of 



,1 -2^1 + Z ') 



(14) 



in terms of the Lyman limit frequency ^ll and x max (n) 
corresponds to the comoving distance between z and 
in) given by: 

[l-(n + I)' 2 ] 



1 + z max (n) = (1 + z)- 



(15) 



(1-n- 2 ) 

If the photon emis sivity, e is homogeneou s, equation [T?] 
can be written as (|Barkana fc Loebll2005f ). 



(-1 , 7 \2 "max 

j a ( Z ) = [ -^ L E 



4~ 



n=1 



* (n) cdz' 



H(z') 



(16) 



3.3. Gas temperature 



Once star formation has got underway a population of 
stellar remnants will be produced capable of generating 
highly energetic X-rays. Several candidate X-ray sources 
exist including X-ray binaries in starburst galaxies, in- 
verse Compton scattering of CMB photons f rom rela- 
tivistic electrons in supernova remnants (SNR) (jOh et alJ 
l200l . and mini-quasars. X-rays may contribute to reion- 
ization, although constraints from the unresolved soft 
X-ray background suggest that t his is not the dom i- 
nant source of ionizing photons (jDiikstra et al.l 12004) . 
More importantly, as X-rays ionize hydrogen they de- 
posit much of their energy as heat. This X-ray preheating 
can easily be sufficient to heat the IGM above the tem- 
perature of the CMB. Although hard X-rays have a mean 
free path comparable to the Hubble size, most of the 
heating turns out to be done by soft X-rays (E < 2 keV) 
which can produ ce significantly inhomogeneou s heating 
at high redshifts (jPritchard fc Furlanettoll2007f) . 

We calcu l ate th e X-ray heating following the model of 
IFurlanetto! (|2006l ). We model X -ray sources with a spec- 
tral distribution function (number of photons per unit 
comoving volume per unit time and frequency) that is a 
power-law with index as- 



hv \i/ 



— as — 1 



(17) 



and the pivot energy is hv§ = 1 keV. We assume emis- 
sion within the band 0.1 keV to 30 keV and set the 
normalization constant Lq by requiring that the inte- 
grated power density P = J dvL (v / v )~ as is 3.4 x 
10 40 /x ergs -1 Mpc -3 when integrated between 0.1 keV 
and 30 keV, but with a highly uncertain constant factor 
fx- This normalization is chosen so that, with fx = l, 
the total X-ray luminosity per unit SFRD (star forma- 
tion rate density) is consistent with that observed in 
starb urst galaxies in the present epoch (see IFurlanetto! 
l2006l for further details). Typical values are as = 1.5 
for starbursts, as = 1.0 for supcrnovac remnants, and 
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as = 0.5 for the X-ray background generated by mini- 
quasars. These spectral indices are consistent with mea- 
sured spectra of known X-ray sources, though it is largely 
uncertain whether low-redshift sources can be fully con- 
sidered as a representation of source spectra at high red- 
shifts. 

We link the total X-ray emissivity (number of photons 
per SFRD per unit comoving volume per unit frequency 
and time) to the star formation rate by 

. / \ ~ f \ ( SFRD(x,z) \ 
ex(x ; .,,)=e x (,)^ MQyr _ lMpe _3 j. (18) 

The X-ray number flux per unit frequency is then 



J A '(x, z, v) 



,0- + z) 



(19) 

where again |x'| is the comoving distance between z and 
z', and v' is the emission frequency at z' corresponding 
to an X-ray frequency v at z 

../ ..(!+*') 



(1 + *) 

The optical depth is given by 
r{z, v, x, x') = 



(20) 



all [n m o m ( V ")+ (21) 
+«HolCTHcl(^") +n H eII^HeIl(^")] ) 

where the integral is along the photon path in proper 
units between emission (x + x') at redshift z' and re- 
ception (x) at redshift z, while v" is the frequency cor- 
responding to v at the redshift along the photon path. 
The c ross-sections for ioni zation are calculated using the 
fits of iVerner et all rflQQfih . 

Finally, the total rate of energy deposition per unit 
volume is 



ex(x, z) 



dv<Ji{y)Jx{x-i z, v)(hv — hv l th ) 



(22) 

where i =HI, Hel, Hell, rij is the number density and 
hv t \ x = E t h is the threshold energy for ionization. To get 
the total heating rate, we multiply this by the fraction of 
energy converte d into heat .fheat, obtained us ing the fit- 
ting formulae of IShull fc van Steenberd (|l985h . We then 
evolve the gas temperature using 

AT K _ 2T K An | 2e x /heat 
At 3n At 3fcsn ' 

where n is the proper number density of all particles. In 
order to evolve this equation through our simulation box, 
we set the initial condition for the gas temperature at 
z = 24.9 to be T K = 14.43 K. The latter temperature is 
derived by assuming the adiabatic cooling of the gas since 
recombination for our fiducial cosmological model. In 
setting this initial condition we also assume that the gas 
is homogeneously cooled to this temperature at z=24.9. 
Note that this is the same initial temperature as in our 
constant-temperature case where we ignore fluctuations 
in X-ray heating of gas, among others. 

As /heat depends on the free electron fraction in the 
IGM, x e , we must also evolve x e using 



Ax e 
At 



ex 



nE, 



th 



(24) 



where /i on is the fraction of energy converted into ion- 
izations which also depends on x e (note that ex also de- 
pends on x e through rii in equation I22p. Note that this 
term is quite important since otherwise we get x e ~ 
and /heat ~ in the IGM from UV ionization while 
/heat ~ 0.01-0.1 whenx e ~ 10- 8 -10~ 4 instead of 0. We 
neglect primary ionizations from X-rays since secondary 
ionizations dominate at large radii from halos and UV 
ionizations dominate at small radii. In making the above 

calcu lation, since /heat oc x\ (|Shull fc van Steenberd 
Il985f ). we approximate it by only considering the ion- 
ization of hydrogen and by setting E t h = 13.6 eV. Both 
recombinations and corrections from helium are typically 
not important for calculating x e , which stays small over 
the redshift range of interest. Figure [2] shows the evolu- 
tion of x e in the IGM for the simulation, which remains 
small at all times. Note that x e is defined in the neutral 
IGM outside of fully ionized HII regions. 

4. THE 21-CM SIGNAL: SIMULATIONS 

In order to obtain the 21-cm brightness temperature 
from simulations, we basically need to apply equation 
[51 Both /hi and pb are already obtained by the simula- 
tion as these properties are calculated with the evolution 
of the dark matter distribution. On the other hand, the 
spin temperature, is calculated at the post-processing 
stage since the radiative transfer calculations of the ini- 
tial run do not take into account the Lya photons and 
X-ray heating. Usually it is assumed that T$ >> T 7 (e.g. 
the number density of hydrogen atoms in the triplet level 
is saturated) so that one does not need to worry about 
the spin temperatur e contribution to the 21-cm signal 
(|Mellema et al.ll2006l ). However, this assumption should 
only be safe for z < 10 so that one needs to consider 
the evolution of the spin temperature for a proper treat- 
ment of the 21-cm signal at the higher redshifts provided 
by this simulation. Moreover, fluctuations in the Lya 
coupling and X-ray heating, may be important at high 
redshifts dBarkana fc Loeb 2005; Prit chard fc Furlanettol 
l2007t ISemelin et al.ll2007t ). Therefore, we present here a 
full calculation of these fluctuations on the high redshift 
signal probed by the simulation. 

4.1. Coupling parameters 

Calculation of the collisional coupling parameter from 
equation [TU] is straightforward. In this case we can easily 
include perturbations due to fluctuations in tt-hi and n e . 
In order to obtain the radiative coupling parameter, y a 
we need to determine the comoving photon emissivity 
directly from the simulation using: 

e Q (x, u, z) = SFRD{ X , z)e b {v) , (25) 

where SFRD(x,z) is again the star formation rate den- 
sity from the simulation (in terms of the number of 
baryons in stars per comoving volume and proper time) 
and tbiy) is the spectral distribution function of the 
sources (defined as the number of photons per unit 
frequency emitted at v per baryon in stars). Note 
that we are assuming that stars dominate over mini- 
quasars for the radiative coupling. We consider sepa- 
rately the spectral dist ribution function fro m Pop II stars 
(iLeitherer et al.lll999h and Pop III stars (jBromm et alJ 
120011 but see also lBarkana fc Loebll2004h . We then ap- 
ply directly equation [14] to the simulation boxes. We can 
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speed up this calculation by noting that the integral can 
be written in terms of a convolution between the SFRD 
and a specified kernel. 

Figure [5l left panel, shows the power spectrum of the 
J a fluctuations from the simulation. The power spec- 
trum is dominated by large scale fluctuations at low red- 
shifts, since most of the photons propagated through the 
entire box and fluctuations represents a small percentage 
of the overall Lya photon flux. At large redshift (around 
z ~ 20), most of the Lya photons have not had time to 
propagate very far from halos and cluster around halos 
in 0.2-0.3 hr 1 Mpc bubbles. However, these bubbles are 
highly clustered and have large bias factors related to the 
dark matter density field. The strong clustering of these 
small bubbles around first-light sources still leave a Lya 
intensity power spectrum that is strongly clustered even 
at high redshifts. Our measurement directly by post- 
processing the star-formation rate in our numerical sim- 
ulation to extract properties of gas physics and 21-cm 
brightness temperature shows that there is no Poisson or 
shot-noise component in the Lya intensity fluctuations. 
Analytical models that motivate detections of first-light 
galaxies from the 21-cm background have suggested a 
large s hot-noise for the Lya in tensity background fluctu- 
ations (jBarkana fc Loebl |2005[), but we do not see such a 
component in our simulations to the extent that we can 
separate the Lya intensity power spectrum at redshifts 
z ~ 20. While we cannot comment on the appearance of 
a shot-noise even at higher redshifts or smaller scales, we 
hope to return to the general issue of detecting signatures 
of the highest redshifts galaxies in the 21-cm background 
in an upcoming paper by fully taking into account the 
instrumental systematics and foregrounds. 

Finally the Lya coupling is determined through equa- 
tion [11] (there will be extra fluctuations due to the S a 
dependence on Tk and the matter density distribution). 
Analysis of the simulation shows that on average y a > 1 
for z < 17 signaling the approach of the spin tempera- 
ture to the gas temperature. Moreover, the Lya coupling 
dominates over collisions up to z ~ 22 when y c ~ 10 -2 . 

4.2. Gas temperature 

As described in Section [2~T1 the gas temperature pro- 
vided by the simulation is essentially due to the virial 
temperature, given by the velocity dispersion of each par- 
ticle. This effectively sets the gas temperature to T > 10 4 
K within the virial radius of halos. However, most of the 
heating is restricted to the high density, ionized regions, 
while the neutral IGM continues to cool adiabatically. 
This means that most of the 21-cm signal would be seen 
in absorption even at the low redshifts when reionization 
is well underway 

X-ray heating on the other hand should heat the neu- 
tral IGM above the CMB temperature fairly easily and 
needs to be taken into account for a proper treatment 
of the 21-cm signal. We need therefore to include X-ray 
heating as part of our post-processing to predict the spin 
temperature of gas and to calculate the brightness tem- 
perature of the 21-cm signal. Due to the clear separation 
between the ionized and neutral regions, we can consider 
the evolution of both heating mechanisms separately. 

In order to calculate heating due to X-rays, we basi- 
cally follow the procedure outlined in section 13.31 using 
the star formation rate density provided by the simula- 



tion. In order to use the same Fourier transform tech- 
nique to speed up the temperature fluctuation calcula- 
tion, we assumed the density of our species to be spatially 
constant and equal to the average in the computation of 
the optical depth. This approximation will reduce the 
optical depth around the halo, but as one get further 
away our optical depth should converge to the real one. 
We therefore probably underestimate the temperature 
in the halos and overestimate in the neighboring regions 
outside. For the spectral distribution function we as- 
sume as = 1.5 (starbursts). Figure[5l right panel, shows 
the dimensionless power spectrum of the gas tempera- 
ture due to X-rays from the simulation. Note that these 
fluctuations are only relevant to the 21-cm signal as long 
as the gas temperature is not much higher than the CMB 
temperature. Thermal histories are also plotted in Fig- 
ure [HI for fx = 0.1, 1.0, and 10.0. These indicate that X- 
ray preheating can indeed heat the gas above the CMB 
temperature at the redshift range important for 21-cm 
observations, justifying the assumption that 3> Tcmb 
at redshifts z < 10. Clearly, though, there is considerable 
uncertainty in what the exact thermal history is expected 
to be. Hereafter, when we calculate the 21-cm brightness 
temperature and its anisotropy, we will assume fx = 1-0. 
The analytical result we present here is also matched to 
the same X-ray intensity. 

To highlight the reason why we consider X-ray heating 
as an inhomogeneous process, in Fig. 7, we plot the vol- 
ume filling factor of our simulation calculated by taking 
the ratio X^(4/37rA 3 )/VJji mu i where A is the mean free 
path around each X-ray source identified in the simula- 
tion with volume V^imui- While photons with energies 
around 100 eV that are primarily responsible for heat- 
ing propagate rapidly and fill the box by z ~ 10, the 
volume filling factor is below 1 around z ~ 15. This 
demonstrates that fluctuations in the heating of the gas 
by X-rays will be important around these redshifts. 

4.3. Brightness temperature 

Finally, implementing equation [51 we can calculate the 
21-cm brightness temperature for the simulation. Figure 
[H] shows the average evolution of this brightness temper- 
ature, together with the CMB temperature, gas temper- 
ature, and the spin temperature. Note that, as already 
pointed out, the spin temperature decouples from the 
CMB at z < 17 when y a > 1. The average 21-cm sig- 
nal is clearly non-negligible at high redshifts. Moreover 
the transition from emission to absorption depends on 
whether X-ray fluctuations are included. This is because 
regions with very low ionization fraction are still cold in 
spin temperature at low redshifts relative to the CMB 
and are therefore seen in absorption. 

Figure [H shows the evolution of the 21-cm signal with 
redshift for a set of k values. Again we can see that 
the signal strength actually increases for z > 12. The 
evolution of the 21-cm brightness fluctuations resemble 
strongly the average 21-cm brightness temperature. The 
fluctuations increase up to z = 15 as the coupling of 
the spin temperature to the gas temperature increases, 
then diminish during the absorption-emission transition 
down to z = 11, increasing or reaching a plateau in the 
Ts >> Tcm b i regime and then falling again at z = 8 — 9 
when the ionization fraction reaches very high values and 
the universe fully reionizes. 



11 




z=7.40 
z=10.0 
z=15.2 
z=20.6 



k [h/Mpc] 



100 
10 

1 

0.1 
' 0.01 
10~ 3 
10-" 
10~ 5 
10~ 6 



k [h/Mpc] 



z=7.40 
z=10.0 
z=15.2 
z=20.6 



Fig. 5. — Left: Power spectra of the Lya flux at several redshifts divided by the square of the mean flux. Right: power spectra of the 
gas temperature due to x-rays divided by the average squared at each redshift. 
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Fig. 6. — Evolution of gas temperature due to X-rays for several 
fx (normalization of the X-ray luminosity). 

The evolution of the power spectrum on large 
scales is qualitat i vely similar to that calculated in 
iPritchard fc Loebl (|2008f ) showing three peaks resulting 
from the periods where ionization, temperature, and Ly- 
alpha fluctuations respectively come to dominate. We 
note that the dip at Xi ~ 0.3 occurs because on large 
scales the gas temperature is quite close to the CMB 
temperature at this redshift. 

Figure [10] shows the brightness temperature for the 
same maps as in figure [U including fluctuations in the 
matter density, ionization fraction, Lya coupling and X- 
ray heating. Note that some of the features from the 
rcionization maps are somewhat smeared in the bright- 
ness temperature maps due to the convolution with the 
density. At z ~ 15 perturbations in Lya and X-ray heat- 
ing are already non-negligible and need to be properly 
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Fig. 7. — Volume filling factor of the X-ray radiation in the 
simulation as a function of redshift for several different values of 
the photon energy. While a photon with an energy such as 100 cV, 
that dominates X-ray heating, fills up the whole volume by z ~ 10 
(ie. mean free path is larger than the simulation box length), at 
z ~ 15, the volume filling factor remains below 1 and X-ray heating 
at the onset of rcionization is inhomogencous. 

taken into account when making 21-cm predictions. 

Figure [IT] shows the 21-cm signal at z = 20, 15 and 
z = 10 from top to bottom respectively, where we sepa- 
rately include fluctuations in following quantities: matter 
density plus ionization fraction (left panel), matter den- 
sity plus ionization fraction plus Lya coupling (middle 
panel), and matter density plus ionization fraction plus 
X-ray heating (right panel). While at z ~ 20 fluctua- 
tions in the Lya intensity field are important, at z ~ 15, 
features in the brightness temperature are sensitive to 
whether we model X-ray heating as an inhomogeneous 
source or whether we consider X-ray heating, wrongly, as 
uniform. This is the redshift in our simulation at which 
, due to X-ray heating, regions are beginning to be first 
detected in emission of the 21-cm brightness temperature 
instead of absorption as at higher redshifts. 

In terms of the inhomogeneous sources of the 21-cm 
brightness temperature, fluctuations in Lya background 
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Fig. 8. — Temperature of the CMB (black solid line), gas (red 
dashed line) and spin (blue dotted line) as a function of 2, in our 
simulation, in which we have included x-ray heating, collisional 
and radiative coupling. The lower, solid black line shows the av- 
erage brightness temperature when all fluctuations are taken into 
account, while the dashed one shows the absolute value (since it 
is negative at these redshifts). Lower, red (grey) line shows the 
brightness temperature when only fluctuations in the matter den- 
sity and ionization fraction are included. 
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Fig. 9. — Evolution of the brightness temperature power spec- 
trum with redshift. 

makes the biggest difference at z ~ 20 as can be seen 
from the bottom middle panel. Here, the locations with 
largest absorption signal in the 21-cm background is asso- 
ciated with first-light sources that have just formed first 
stars. A comparison of the middle panel of Fig. [IT] at 
z = 20 and the right panel at z = 15 shows that regions 
that are first dominating in emission due to X-ray heat- 
ing at z ~ 15 are mostly the same regions that were first 
brighter in absorption at z ~ 20 due to Lya coupling. At 
lower redshifts, however, during partial reionization, fluc- 
tuations in the brightness temperature are dominated by 
fluctuations in the ionization fraction modulated by the 
density field inhomogeneities. Thus, at lower redshifts 
z < 10 that will be targeted by the first-generation 21-cm 



interferometers, one can mostly ignore the effects associ- 
ated with inhomogeneous X-ray heating or anisotropics 
in the Lya coupling. This becomes clear when we com- 
pare the predictions related to the 21-cm brightness tem- 
perature anisotropy power spectrum with those related 
to analytical models in the next Section. 

5. OVERALL COMPARISON OF THE 21-CM SIGNAL: 
SIMULATION AND ANALYTICAL MODEL 

Here we compare the power spectrum of the 21-cm sig- 
nal from the simulation to a fast analytical power spec- 
trum generator. The calculation of the analytical power 
spectra will basically f o llow t he procedure outlined in 
iPritchard fe Furlanettol (|2007l ). using the model for the 
ionization fra ction discussed in sectionl2.2| a nd we refer 
the reader to IPritchard fc Furlanettol (|2007f ) for further 
details. 

When considering fluctuations in the 21-cm brightness 
temperature arising from variations in ionization and 
density, the correlation function of the brightness tem- 
perature can be expressed as 



Z TbTb ee ({ST b - ST b ){ST b - Sf b )) 
= T c [fm€ss + Caw (1 + €ss) 



(26) 
,-)] 



and we use the procedure described in section 
to calculate £$8, £x t Xi an d £,x.iS- Note that we as- 
sume gaussianity in this calculation but take into ac- 
count the Gaussian terms from the 4-point function 
(see the appendix for a discussion of the effect of the 
non-Gaussian terms in equation I26jl . However, varia- 
tions in the ionization fraction and density are not the 
only sources of 21-cm brightness temperature fluctua- 
tions and we include anisotropics in X-ray heating and 
Lya coupling separately (by following the procedure in 
IPritchard fc Furlanett"ofl2007f l so we can study how pre- 
dictions change with whether one considers heating to be 
uniform, for example. We also follow the same procedure 
with the analytical calculation by adding various sources 
in step for easy comparison with numerical simulations. 

In order to calculate the gas temperature (due to heat- 
ing by X-rays) , we need to obtain the analytical star for- 
mation rate and use it in equation [TU We model the star 
formation rate as tracking the collapse of matter, so that 
we may write the star formation rate per (comoving) unit 
volume 



SFRD = p°(z)/,-/ coll (z). 



(27) 



where is the cosmic mean baryon density today. This 
formalism is appropriate for z > 10, as at later times 
star formation as a result of mergers becomes important. 
For the analytical calculation, we do not distinguish be- 
tween Pop II and Pop III stars and so use a value of 
/* = 0.1, appropriate for Pop II stars, which dominate 
star formation at lower redshifts. While these param- 
eters have not been fitted to the simulation data, the 
star-formation rates from theory and simulation agree 
quite well. Fourier transforming the correlation function, 
£,T b T b , in equation |2"61 yields the desired power spectra. 
By first generating the correlation functions and then 
Fourier transforming we avoid having to consider the 
power spectrum convolution for the £ XiXi €s6 an d £,xi8^xi8 
terms. 
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Fig. 10. — Maps of the 21-cm brightness temperature from the simulation, including all fluctuations, at redshifts z 
as labeled, corresponding to x t = 0.0002, 0.03, 0.35, 0.84. 
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In figure [TS] we show a comparison of the brightness 
temperature power spectrum between the analytical cal- 
culation and numerical simulations. Previously we have 
commented on the comparison between the analytical 
model prediction for the ionization fraction power spec- 
trum (relative to the density field power spectrum) and 
the numerical simulation for the same quantity (Fig. 3 
and Section 2.2). We also expect the 21-cm prediction 
which uses only the ionization fluctuations and the den- 
sity field as inhomogeneous sources to agree with the 
analytical calculation to the same extent that the two 
agreed previously. Thus, any differences beyond that of 
Fig. 3 in Fig. 1121 between simulations and the analytical 
model are due to differences in the two prescriptions re- 
lated to fluctuations in the Lya intensity field and X-ray 
heating. 

At low redshifts, we note that the two predictions agree 
very well to the extent that we can trust our simula- 
tions. At z ~ 7, the difference at k ~ 0.1 h Mpc -1 is 
most likely to be that of the finite volume of the box 
and the cut-off associated with large-scale bubbles that 



have grown beyond the volume of the simulations, since 
we find less power in the simulation compared to that 
of the analytical model (see figure [3|. At z ~ 10 the 
agreement between the simulation and analytical curves 
corresponding to density plus ionization fraction fluctu- 
ations (dotted lines) is very good. The difference in the 
solid lines is because x-ray perturbations are actually im- 
portant at these redshifts. When we include fluctuations 
in x-rays, the very large, neutral scales (k < 0.4) are still 
"cold" at z = 10 (close to the CMB temperature) so we 
get a reduction in the 21-cm perturbations (this can also 
be seen in Figure \§§ . The differences between the two 
prescriptions become more obvious at z ~ 15 and higher. 
At these high redshifts, as we have already noted, differ- 
ent sources of fluctuations make different contributions 
to the 21-cm brightness temperature fluctuations. 

At z ~ 15, analytical calculations suggest a charac- 
teristic scale for 21-cm brightness temperature fluctua- 
tions with a shot-noise like power-spectrum at k > a 
few h Mpc . This 21-cm power spectrum is dominated 
by fluctuations in the X-ray heating background as can 
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Fig. 11. — Maps of the 21-cm brightness temperature. Top: z = 10. Middle: z = 15. Bottom: z = 20. Left: fluctuations in density 
and ionization fraction with homogeneous X-ray heating and Lya radiation field. Center: Fluctuations in the Lya background added, but 
homogeneous X-ray heating. Right: fluctuations in X-ray heating, but uniform Lya radiation field. 



be established from a comparison between the dotted 
and dashed lines related to the analytical calculation at 
z ~ 15. At this redshift, the analytic model shows a con- 
siderable deficit of power on small scales when compared 
to the simulation and inhomogeneous X-ray heating is 
included. This appears to be a result of the analytic 
model assuming a tight cross-correlation between tem- 
perature and density fluctuations on all scales. Since the 
Pts term contributes with negative sign during the ab- 
sorption epoch, P21 is reduced. It appears however that 
on small scales the cross-correlation between density and 
temperature is small. On large scales, the negative sign 
of the density-Lya cross-correlation can be seen where 
the full calculation lies below than that where Lya fluc- 



tuations are ignored. If we set Pts to zero the analytic 
calculation gives much better agreement with the simula- 
tion on small scales. Clearly, im provements to the simple 
model of X-ray heating given in iPritchard fc Furlanettol 
(|2007| ) are necessary to resolve this problem. Note also 
that the dashed line is slightly higher than the total 
contribution shown by the solid line due to the anti- 
correlation between X-rays and Lya fluctuations. 

At z ~ 20, fluctuations in the X-ray heating are not the 
dominant source of 21-cm brightness temperature fluc- 
tuations, but rather it is a combination of density per- 
turbations and Lya coupling. While there are general 
differences at z > 10 between the analytical model and 
our numerical simulation, we note that these still lead to 
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predictions that agree within a factor of a few, while at 
z < 10, the agreement is better than 30% and at z ~ 7 it 
is significantly better over the wavenumber range where 
we can make a comparison. Given that simple fitting 
formulae can be written down to model the ionization 
fraction power spectrum, which we outlined in Section 
2.2, it is likely that one can quickly calculate the 21- 
cm power spectrum at low redshifts as a function of the 
density field power spectrum. For reasons that we have 
previously discussed to enable quick exploration of the 
parameter space of 21-cm experiments, we suggest that 
analytical calculations, complemented by well-calibrated 
fitting formulae, can be trusted. For parameter estima- 
tion with data from upcoming interferometers, it may be 
necessary to improve beyond the fitting functions, how- 
ever. 

We note that comparisons such as the one we have 
performed between an analytical model and a numeri- 
cal simulation of reionization, which was post-processed 
to extract properties of the 21-cm brightness tempera- 
ture, have been performed in the literature and lead- 
ing to conclusions similar to ours that analytical calcu- 
lations are adequate for parameter predictions and mea- 



surements with 21-cm data (jZahn et al.ll2007[ K Our work 
differs from these previous comparisons in that, for the 
first time to the extent we are aware of, we extend the 
comparisons to redshifts beyond 8 and comment on the 
agreement even out to z ~ 20. To do this, we are forced 
to model fluctuations in X-ray heating and Lya intensity 
field since as we have shown these inhomogeneous sources 
make significant contributions to 21-cm brightness tem- 
perature fluctuations at z > 10. At low redshifts, even 
after including fluctuations in Lya coupling and X-ray 
heating, we find that the fluctuations in the 21-cm bright- 
ness temperature are dominated by fluctuations in the 
ionization fraction and the density field. Thus, previous 
analytical and numerical simulation comparisons, which 
only considered physics on ionization bubbles and their 
distribution in detail, remain valid. 

6. SUMMARY AND CONCLUSIONS 

Using a new large volume, high resolution simulation of 
cosmic reionization based on a hybrid code for N-body 
dark matter and radiative transfer of ionizing photons 
through an adaptive algorithm, we have measured sev- 
eral properties of the reionization process. We have fo- 
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cuscd our discussion on the low- frequency 21-cm signal 
associated with the neutral hydrogen distribution which 
is now pursued by a variety of interferometers as a probe 
of the rcionization history of the universe. 

In this paper we have studied the extent to which 
statistical results from an analytical model arc consis- 
tent with results extracted from the simulation. Note 
that first principle simulations cannot be considered for 
baryon physics and star formation so that the brightness 
temperature of the simulation 21-cm signal is derived by 
post-processing the simulation with certain results based 
on analytical prescriptions of the rcionization process. 

In detail, we have compared the spatial clustering of 
the neutral gas fraction, ionization fraction, and the as- 
sociated 21-cm signal from the neutral hydrogen distri- 
bution. Our study extends to high rcdshifts where the 
contribution from spin temperature is non-negligible and 
we take into account the heating of the gas by X-rays and 
the effect of Lya and inhomogencous collisional coupling 
when calculating the 21-cm radio signal. We find very 
good agreement between simulations and an analytical 
model at low redshifts, though there are, non-negligible 
differences at higher redshifts (z > 10) arising from dif- 
ferences related to X-ray heating and fluctuations in the 
Lya coupling. At the redshift range that will be probed 
by the first-generation 21-cm experiments (z < 9), we 
find that simple analytical models coupled with fitting 
functions associated with the ionization fraction power 



spectrum can easily reproduce the results from the nu- 
merical simulation. Therefore, at these redshifts, we find 
that there are no remaining issues with using an analyt- 
ical model to explore the parameter space relevant for 
future 21-cm surveys, when using estimators based on 
the power spectrum alone. At higher rcdshifts, detailed 
comparisons, especially with regard to gas temperature 
and Lya coupling, may be desirable in order to improve 
current analytical models. 
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APPENDIX 

When non-Gaussian terms become important, cquation[2H]is no longer valid and we need to take in to account the fu ll 
4-point and 3-point function in the power spectrum calculation of the brightness temperature (see lLidz et aLlfeOOTf ). 
The full power spectrum is then 

P2i(k) = Tl [f^PsAk) + Px u xM - VmPxM + 2P XiS , Xi (k) - 2f m P XiS j(k) + P XiS , Xi5 (k)] , (1) 

where P a {, is just the power spectrum between the quantity a and b. 

To see the difference, we show in figure Q2] the power spectrum of the brightness temperature obtained directly from 
the simulation (only considering fluctuations in the density and ionization fraction) , compared to the one obtained by 
just using the first three terms in the equation above (the "low order" terms). Note that all the power spectra used are 
obtained directly from the simulation. We also plot the contribution from the higher order terms (second line in the 
equation) plus the result of considering the full expression, which, as expected, is similar to the actual 21-cm power 
spectrum measured from the simulation. 

We can see that, only at higher rcdshifts can the contribution of the "higher order" terms be safely neglected. In 
addition to accounting for some of the differences in the power spectrum, non-Gaussianity can also be an important 
source of information, specially at the height of the reionization process (i.e., about x~i ~ 0.5), when the ionization 
structure becomes quite non-Gaussian. 
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